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Abstract 

A novel model of biological organisms is advanced, treating an organism as a self-consistent 
system subject to a pathogen flux. The principal novelty of the model is that it describes not 
some parts, but a biological organism as a whole. The organism is modeled by a five-dimensional 
dynamical system. The organism homeostasis is described by the evolution equations for five 
interacting components: healthy cells, ill cells, innate immune cells, specific immune cells, and 
pathogens. The stability analysis demonstrates that, in a wide domain of the parameter space, the 
system exhibits robust structural stability. There always exist four stable stationary solutions char- 
acterizing four qualitatively differing states of the organism: alive state, boundary state, critical 
state, and dead state. 
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1 Introduction 



Biological organisms are among the most complex systems, requiring for their mathematical descrip- 
tion rather elaborate equations [1-3]. Such equations are, as a rule, nonlinear, because of which 
they can possess different solutions, depending on the system parameters. The solutions of nonlin- 
ear equations are known to be very sensitive to the values of the control parameters. Sometimes, 
a slight variation of a parameter results in a discontinuous change of the system behavior, which is 
manifested in the qualitative change of its phase portrait. Such abrupt transformations of the phase 
diagram demonstrate the structural instability of the related dynamical system and are also called 
catastrophes [4]. In many cases, these catastrophes describe a real morphogenesis occurring in com- 
plex systems under the variation of the parameters. However, structural instability may also happen 
as an artifact, simply because a real complex system has not been correctly modelled by a dynamical 
system. It is therefore extremely important to carefully take into account the basic features of the 
real system, when modelling it by mathematical equations. Even a small term can essentially influ- 
ence the behavior of a nonlinear dynamical system [5]. This especially concerns the terms reflecting 
fundamental symmetries of the complex system. 

In the present paper, we consider such a very complicated biological system as an organism con- 
sisting of five components, healthy cells, ill cells, two types of immune cells, innate and specific, and 
pathogens. Our aim is twofold. First, we formulate the dynamical system describing the homeostasis 
of an organism as a whole, paying attention to the necessity that the action-counteraction exchange 
symmetry be preserved. The principal point in our description is the treatment not of some parts of an 
organism, but the study of the latter as a self-consistent system. Second, we perform the stability anal- 
ysis and demonstrate that the system exhibits a remarkable structural stability: in a very wide range 
of parameters, there always exist only four stable stationary solutions characterizing four organism 
states: alive state, boundary state, critical state, and dead state. Such a structural stability found in the 
proposed dynamical system provides an important validation step for the model [6], as it captures one 
of the most important characteristics of biological organisms, that of adapting robustly to different 
conditions. 

2 Construction of evolution equations 

Let us consider biological species, enumerated with the index i = 1,2,..., the number of agents of 
the z-th type being iVj. Evolution equations can be represented either by differential or difference 
equations. Keeping in mind a very large number of interacting agents and the interaction times short 
as compared to the observation time, we employ here a picture with continuous time. The general 
structure of the dynamical system, describing the evolution of the species, can be represented as a set 
of differential equations 



where t is time, Ri is a life rate, Rij is an interaction intensity, and Fj is an influx. The quantities Ri 
and R,^ are treated as parameters. The influx F describes an input that is external with respect to the 
z-type species. That is, Fj may include a flux that is external for the organism as a whole and also it 
may describe a flux from other species. Generally, F can be a function of Nj, with j ^ i, such that 
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To be self-consistent, the set of equations © has to satisfy two major requirements. First of all, 
the various processes included in the consideration must be of the same order of nonlinearity. This 
means that, since the interaction term in Eq. (QQ) is of second order, the fluxes Fi have also to be not 
higher than of second order with respect to Nj, with j ^ i. 

Another important requirement is the existence of the action-counteraction symmetry. This means 
that, if in Eqs. (OQ) there occurs a term RijNiNj then, these equations have to contain the exchange term 
RjiNjNi. That is, if there is an action, there should exist a counteraction, which can be schematically 
formulated as the symmetry 

RijNiNj <-> RjiNjNi . (3) 

Note that Rij ^ Rji in general. 

To specify a biological system, we consider an organism consisting of Ni healthy cells, N 2 ill 
cells, N 3 innate immune cells, N 4 specific immune cells, and N 5 pathogens. The organism home- 
ostasis is characterized by the processes and reactions between the cells and pathogens. A detailed 
medical description of these processes and reactions is given in Ref. [7]. Here, we shall list them 
only in brief, specifying the parameters Rj and R^, with taking account of the typical signs of these 
parameters. 

(i) Healthy cells are characterized by a natural reproduction rate R\ = A\. Since the body volume 
is limited, the carrying capacity limitation is governed by Rn = —An. Ill cells do not interact directly 
with healthy cells, which implies that R\ 2 = 0. Healthy cells can be occasionally attacked by innate 
as well as specific immune cells, hence R\ 3 = — A 13 and R u = —Au, which causes autoimmune 
diseases. Pathogens infect healthy cells, thus R 15 = — v4 15 . We treat an organism as a self-organizing 
system, which means that healthy cells can be reproduced inside the organism, but are not supplied 
from exterior, that is, F\ = 0. 

(ii) /// cells have either a natural death rate R 2 = — A 2 , when A 2 > 0, or can exhibit an unnatural 
proliferation, when A 2 < 0, which happens under certain diseases, such as cancer. Healthy cells do 
not directly interact with ill cells, so R 2 \ = 0. The carrying capacity limitation for ill cells, while it 
may exist, is considered to be much larger than that for healthy cells, which allows us to set R 22 = 0. 
Ill cells are killed and eliminated by immune cells, hence R 23 = —A 23 and R 24 = —A 24 . The 
degradation of ill cells is increased under the influence of pathogens as the latter catalyze the immune 
system, hence R 2 $ = —A 25 . The number of ill cells rises as a result of pathogens infecting healthy 
cells, which gives F 2 = A 51 N 5 Ni. 

(iii) Innate immune cells die by apoptosis, with a rate R 3 = —A 3 . They can be promoted by 
healthy cells, so R 3 i = A 31 . And they are activated by ill cells, R 32 = A 32 . The carrying capacity 
of immune cells is much larger than that of healthy cells, which makes it admissible to set R 33 = 0. 
Innate immune cells are activated by specific immune cells, hence i? 34 = A 34 . And they are activated 
by pathogens, R 3 $ = A 35 . There is no external flux from outside of the organism, that is F 3 = 0. 

(iv) Specific immune cells also have a finite lifetime, characterized by an apoptosis rate R 4 = — A 4 . 
They are promoted by healthy cells, 7? 4 i = A41, and are activated by ill cells, i? 42 = A i2 . Innate 
immune cells inhibit an excessive amount of specific immune cells, R$ 3 = —^43. Similarly to innate 
cells, for specific immune cells, the carrying capacity limitation can be ignored, so that Ru = 0. 
Pathogens activate specific immune cells, hence _R 45 = A i5 . And there is no external flux, F 4 = 0. 

(v) Pathogens are characterized by a natural decay rate R 5 = —A 5 . Their number does not depend 
on the number of healthy cells, R 5 i = 0. Pathogens proliferate by the lysis of ill cells, R 52 = A 52 . 
They are killed and eliminated by innate, as well as specific immune cells, hence R 53 = — A 53 and 
i? 54 = — t4 54 . The number of pathogens can be an order of or several orders larger than that of healthy 



cells. Therefore, there is practically no carrying capacity limitation for them, that is, it is safe to set 
i?55 = 0. Contrary to all other organism cells, pathogens are supplied from the external surrounding, 
therefore F 5 = F is not zero, but is to be treated as a parameter characterizing the environment in 
which the organism lives. In this letter, we only focus on a constant pathogen flux. 

Taking into consideration the described processes transforms Eq. (OQ) to the set of five evolution 
equations 
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-A 2 N 2 - A 23 N 2 N 3 - A 24 N 2 N 4 - A 25 N 2 N 5 + A 51 N 5 N X , 
-A 3 N 3 + AnNsN! + A 32 N 3 N 2 + A 34 N 3 N 4 + A 35 N 3 N 5 , 
-A 4 N 4 + A 41 N 4 N X + A 42 N 4 N 2 - A 43 N 4 N 3 + A 45 N 4 N 5 , 



-A 5 N 5 + A 52 N 5 N 2 - A 53 N 5 N 3 - A 54 N 5 N 4 + F . (4) 
at 

Using explicitly Eqs. © is not convenient, since the numbers of cells are extremely large. For 
instance, the number of healthy cells is Ni ~ 10 13 , and the number of pathogens can be as large 
as A^ 5 ~ 10 14 or more. Therefore, to be practical, Eqs. © are to be normalized by introducing a 
normalization constant N and defining the cell fractions 

Xi = ^r (z = l,2,3,4,5). (5) 

Note that the cell fractions do not sum up to 1, since N is arbitrary (but is introduced so that we deal 
with quantities of order 1). 

Also, it is necessary to determine a time scale characterizing a typical duration of the homeostasis 
processes. We thus introduce a typical time scale r. Then we define the dimensionless decays rates 

Ok = AiT (6) 

and the dimensionless interaction parameters 

Oij = AijNr . (7) 
In addition, we define the dimensionless pathogen influx 

^^F. (8) 

Making use of the dimensionless quantities, and measuring time in units of r, we reduce Eqs. © 
to the dynamical system 

d t ■ 

^ = fi (<= 1,2,3,4,5), (9) 

with the right-hand sides 

fl = Ot\X\ — CL\\x\ — d\ 3 X\X 3 — CL\ 4 X\X 4 — 0.15X1X5 , 

f 2 — —oi 2 x 2 — a 23 x 2 x 3 — a 24 x 2 x 4 — ci 2 ^x 2 x^ + a^x^Xi , 
h = -a 3 x 3 + a 31 x 3 x 1 + a 32 x 3 x 2 + a 34 x 3 x 4 + a 35 x 3 x 5 , 

f 4 = — OL 4 X 4 -\~ CL 4 iX 4 X\ ~\- Cl 42 X 4 X 2 — d 43 X 4 X 3 -\- G45X4X5 , 

fb = -«5^5 + ^52^5^2 - ^53^5^3 - 054^5^4 + V ■ (1°) 

Equations © and (TTOl) are the basic equations describing the organism homeostasis. 



3 Structural stability analysis 



The dynamical system, given by Eqs. © and (flOl) . contains 25 parameters aj and a^, which appears 
to be an insuperable obstacle for studying its properties. The situation can be simplified by choosing 
appropriate scaling parameters N and r. Thus, for the normalization constant N, we can take the 
capacity number of healthy cells 

N = ^. (ID 
As the temporal scale r, it is reasonable to choose the characteristic time of healthy-cell reproduction 

t=±. (12) 

With these scaling parameters, we have 

«i = 1 , an = 1 . (13) 

These values characterize the typical rates and interactions involved in the homeostasis of the organ- 
ism. 

As is noted above, ill cells can either exhibit a natural decay, when a 2 > 0, or can show a 
pathological proliferation, when a 2 < 0. It is convenient to use the notation 

= , a 2 = 1 - 2/3 . (14) 

The value (3 = corresponds to a 2 = 1 > 0, while (3=1 corresponds to a 2 = — 1 < 0. 
We may assume that the innate and specific immune cells enjoy the same apoptosis rate 

a ee a 3 = a 4 . (15) 

And let us set = 1. We denote the parameters, associated with the interactions between immune 
and healthy cells as 

b = a 13 = a 3 i = a w = a 4 i • (16) 

When 6 = 0, immune cells do not attack healthy cells, hence there are no autoimmune diseases. 
Conversely, for b > 0, autoimmune disorders become possible. We also may assume that all other 
processes, except those related to Eq. ([Tot , are characterized by the intensities being comparable to 
an (equal to 1 according to Eq. (fT3V). That is, we can put 



a 



[a i5 ± b) . (17) 



In this way, we are left with four parameters, a, (3, b, and (p. 

The values of the parameters (3 and b control the occurrence of a chronic pathology or of an 
autoimmune disorder. Varying these parameters, we can reach four limiting cases. 

(i) No chronic pathology and no autoimmune disorder. 

(3 = 0, 6 = 0. (18) 



(ii) No chronic pathology with autoimmune disorder. 



(3 = , 6=1 



(19) 



(iii) Chronic pathology but no autoimmune disorder. 



P=l, 6 = 0. (20) 

(iv) Chronic pathology and autoimmune disorder. 

13 = 1, 6=1. (21) 

We can change the parameters [3 and b in the range between these limiting cases, thus, considering 
a variety of different organisms. For each fixed pair of (3 and b, we vary the parameters a and (p in 
a wide range, finding all admissible fixed points of the dynamical system. In the standard way, we 
accomplish the Lyapunov stability analysis and select the stable fixed points. In order to present the 
results of our analysis in the most illuminating form, it is useful to define the sum of the healthy-cell 
fraction and of the ill-cell fraction, 

x = x\ + x 2 . (22) 
Also, we introduce the summary fraction of all immune cells 

y = x 3 + x 4 . (23) 

The main, and to some extent surprising, conclusion of the stability analysis for the dynamical 
system © is its remarkable structural stability. Varying the parameters (3 and b between the qualita- 
tively different cases (TT8T ) to (I2TI) always results, for any given (3 and b, in four stationary organism 
states characterized by the fixed-point values x* and y*. 

A. Alive state: 

x* > , y* > , (24) 
when there are both self-cells and immune cells. 

B. Boundary state: 

x * > , y* = , (25) 
when there are self-cells, but no immune cells. 

C. Critical state: 

x* = , y* > , (26) 

when only immune cells survive. 

D. Dead state: 

x* = , y* = , (27) 

when there are neither self-cells nor immune cells. 

All these states are always present in the phase diagram on the a — ip plane. The boundaries 
between the states, of course, are different for different (3 and b, but all four states do remain stable. 



The exact domains of stability for each of the states (1241) to (1271) . as well as the values of all 
stationary fractions x*, have been found numerically and, in some cases, analytically. Here, we 
illustrate the results by the phase portraits corresponding to the limiting cases (TT81) to (I2TT) . 

The case (TT81) of no chronic pathology and no autoimmune disorder is represented in the phase 
portrait of Fig 1. The case (fT9l) of no chronic pathology but with autoimmune disorder is shown 
in Fig 2. In Fig. 3, the case (1201) is demonstrated, when there is chronic pathology but there is no 
autoimmune disorder. And Fig. 4 illustrates the case (T2TI) . when there exist both chronic pathology as 
well as autoimmune disorder. As is seen, the boundaries between the states move when varying the 
system parameters, but all four states ([24]) to (l27l) are always present. The case (1201) differs from other 
cases by the presence of a narrow region of bistable states. In this region, there exist large fluctuations 
of the coexisting states [8,9]. The dead state ([27]) is always located in the same part of the phase 
diagram. 

We have also studied the case, when the pathogen influx randomly fluctuates around its mean 
value ([8]). In that case, the cell fractions also fluctuate in time. Then the picture, presented here, 
corresponds to the behavior of guiding centers that are defined by the averaging method [10]. 

Concluding, we have constructed a dynamical system describing the homeostasis of an organism 
as a whole. The organism is composed of five types of components: healthy cells, ill cells, innate 
and specific immune cells, and pathogens. A principal novelty in formulating the model is that we 
consider not some parts of an organism, but treat it as a total self-consistent system subject to the 
influence of pathogens. Another important point in constructing the evolution equations is that we take 
into account the action-counteraction symmetry ©. By varying the system parameters in a very wide 
range, we have shown that the dynamical system © enjoys a remarkable structural stability, always 
exhibiting four stable stationary states. These results remain valid if we include in the dynamical 
system © nonzero values of the carrying capacity limitations for ill cells and immune cells. We have 
also accomplished direct numerical solution of Eqs. ([9]) and ([10] ), confirming the structural stability 
of the dynamical system. 

Our basic aim here has been to present the new model of a biological organism and to accomplish 
a detailed stability analysis making it possible to find numerically the typical phase portraits. An 
important result of the present paper is the remarkable structural stability of the studied biological 
system. The found four stable states characterize all basic states of an organism: alive relatively 
healthy state, boundary ill state, critically ill state, and the dead state. 

We do not overload this paper by discussing various admissible medical interpretations and appli- 
cations. This analysis will be the topic of separate publications. 
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Figure captions 



Fig. 1. Phase portrait on the a — ip plane for the case ([T8T ) of no chronic pathology ((3 = 0) and 
no autoimmune disorder (b = 0), showing the stability domains of the stationary states (1241) to (127T ). 

Fig. 2. Phase portrait for the case (fT9b of no chronic pathology (/3 = 0) but with autoimmune 
disorder (b = 1), demonstrating the stability regions of the stationary states (l24l) to (|27l) . 

Fig. 3. Phase portrait for the case (|20l) . when there exists chronic pathology (J3 = 1) but there is 
no autoimmune disorder (b = 0). The peculiarity here is in the occurrence of bistable states. 

Fig. 4. Phase portrait for the case (I2TT) . when both chronic pathology (J3 = 1) and autoimmune 
disorder (b = 1) are present. The classification of the stationary states is the same as in Eqs. (|24l) to 
(El). 




Figure 1 : Phase portrait on the a — cp plane for the case (fT8l) of no chronic pathology {(3 = 0) and no 
autoimmune disorder (b = 0), showing the stability domains of the stationary states (1241) to (127T ). 




0.0 0.4 0.8 1.2 1.6 2.0 



Figure 2: Phase portrait for the case (fT9l) of no chronic pathology (J5 = 0) but with autoimmune 
disorder (b = 1), demonstrating the stability regions of the stationary states (l24l) to (1271) . 
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Figure 3: Phase portrait for the case (|20l) . when there exists chronic pathology ((3 = 1) but there is no 
autoimmune disorder (b = 0). The peculiarity here is in the occurrence of bistable states. 




Figure 4: Phase portrait for the case (|2D) . when both chronic pathology (J3 = 1) and autoimmune 
disorder (b = 1) are present. The classification of the stationary states is the same as in Eqs. (|24l) to 



